module radiation

  !-----------------------------------------------------------------------------
  ! CAM interface to the legacy 'camrt' radiation code
  !-----------------------------------------------------------------------------

  use shr_kind_mod       , only: r8=>shr_kind_r8, cl=>shr_kind_cl
  use spmd_utils         , only: masterproc
  use ppgrid             , only: pcols, pver, pverp, begchunk, endchunk
  use physics_types      , only: physics_state, physics_ptend
  use phys_grid          , only: get_ncols_p
  use camsrfexch         , only: cam_out_t, cam_in_t    
  use physconst          , only: cpair, cappa
  use time_manager       , only: get_nstep, is_first_restart_step, &
                                 get_curr_calday, get_step_size
  use cam_control_mod    , only: lambm0, obliqr, mvelpp, eccen
  use radae              , only: abstot_3d, absnxt_3d, emstot_3d, initialize_radbuffer, ntoplw
  use scamMod            , only: scm_crm_mode, single_column, have_cld, cldobs, &
                                 have_clwp, clwpobs, have_tg, tground
  use cam_grid_support   , only: cam_grid_write_attr, cam_grid_id,            &
                                 cam_grid_header_info_t, cam_grid_dimensions, &
                                 cam_grid_write_dist_array, cam_grid_read_dist_array
  use cam_history        , only: outfld, hist_fld_active
  use cam_history_support, only: fillvalue
  use cam_pio_utils      , only: cam_pio_def_dim
  use pio,               , only: file_desc_t, var_desc_t, &
                                 pio_double, pio_int, pio_noerr, &
                                 pio_seterrorhandling, pio_bcast_error, &
                                 pio_inq_varid, &
                                 pio_def_var, pio_def_dim, &
                                 pio_put_var, pio_get_var
  use cam_abortutils     , only: endrun
  use error_messages     , only: handle_err
  use perf_mod           , only: t_startf, t_stopf
  use cam_logfile        , only: iulog

  implicit none

  private

  save

  public radiation_readnl
  public radiation_register
  public radiation_nextsw_cday
  public radiation_do
  public radiation_init
  public radiation_define_restart
  public radiation_write_restart
  public radiation_read_restart
  public radiation_tend
  public rad_out_t

  type rad_out_t
    real(r8) solin   (pcols)          ! Solar incident flux
    real(r8) fsntoa  (pcols)          ! Net solar flux at TOA
    real(r8) fsutoa  (pcols)          ! upwelling solar flux at TOA
    real(r8) fsntoac (pcols)          ! Clear sky net solar flux at TOA
    real(r8) fsnirt  (pcols)          ! Near-IR flux absorbed at toa
    real(r8) fsnrtc  (pcols)          ! Clear sky near-IR flux absorbed at toa
    real(r8) fsnirtsq(pcols)          ! Near-IR flux absorbed at toa >= 0.7 microns
    real(r8) fsntc   (pcols)          ! Clear sky total column abs solar flux
    real(r8) fsnsc   (pcols)          ! Clear sky surface abs solar flux
    real(r8) fsdsc   (pcols)          ! Clear sky surface downwelling solar flux
    real(r8) flut    (pcols)          ! Upward flux at top of model
    real(r8) flutc   (pcols)          ! Upward Clear Sky flux at top of model
    real(r8) flntc   (pcols)          ! Clear sky lw flux at model top
    real(r8) flnsc   (pcols)          ! Clear sky lw flux at srf (up-down)
    real(r8) fldsc   (pcols)          ! Clear sky lw flux at srf (down)
    real(r8) flwds   (pcols)          ! Down longwave flux at surface
    real(r8) fsnr    (pcols)
    real(r8) flnr    (pcols)
    real(r8) fsds    (pcols)          ! Surface solar down flux
    real(r8) fln200  (pcols)          ! net longwave flux interpolated to 200 mb
    real(r8) fln200c (pcols)          ! net clearsky longwave flux interpolated to 200 mb
    real(r8) fsn200  (pcols)          ! fns interpolated to 200 mb
    real(r8) fsn200c (pcols)          ! fcns interpolated to 200 mb
    real(r8) sols    (pcols)          ! Solar downward visible direct  to surface
    real(r8) soll    (pcols)          ! Solar downward near infrared direct  to surface
    real(r8) solsd   (pcols)          ! Solar downward visible diffuse to surface
    real(r8) solld   (pcols)          ! Solar downward near infrared diffuse to surface
    real(r8) qrsc    (pcols,pver)     ! clearsky shortwave radiative heating rate
    real(r8) qrlc    (pcols,pver)     ! clearsky longwave  radiative heating rate
    real(r8) fsdtoa  (pcols)          ! Solar input = Flux Solar Downward Top of Atmosphere
    real(r8) swcf    (pcols)          ! shortwave cloud forcing
    real(r8) lwcf    (pcols)          ! longwave cloud forcing

    real(r8) tot_cld_vistau (pcols,pver)  
    real(r8) tot_icld_vistau(pcols,pver) 
    real(r8) liq_icld_vistau(pcols,pver)  ! in-cld liq cloud optical depth (only during day, night = fillvalue)
    real(r8) ice_icld_vistau(pcols,pver)  ! in-cld ice cloud optical depth (only during day, night = fillvalue)
  end type rad_out_t

  character(cl) :: absems_data
  integer :: iradsw      = -1             ! freq. of shortwave radiation calc in time steps (positive) or hours (negative).
  integer :: iradlw      = -1             ! frequency of longwave rad. calc. in time steps (positive) or hours (negative).
  integer :: iradae      = -12            ! frequency of absorp/emis calc in time steps (positive) or hours (negative).
  integer :: irad_always = 0              ! Specifies length of time in timesteps (positive) or hours (negative) SW/LW radiation will be
                                          ! run continuously from the start of an initial or restart run
  logical :: use_rad_dt_cosz = .false.    ! if true use zenith angle averaged over interval between radiation calculations

  integer :: qrs_idx      = 0
  integer :: qrl_idx      = 0 
  integer :: fsds_idx     = 0
  integer :: fsns_idx     = 0
  integer :: fsnt_idx     = 0
  integer :: flns_idx     = 0
  integer :: flnt_idx     = 0
  integer :: cld_idx      = 0
  integer :: rel_idx      = 0
  integer :: rei_idx      = 0
  integer :: cicewp_idx   = -1
  integer :: cliqwp_idx   = -1
  integer :: cldemis_idx  = -1
  integer :: cldtau_idx   = -1
  integer :: nmxrgn_idx   = -1
  integer :: pmxrgn_idx   = -1

  ! Averaging time interval for zenith angle
  real(r8) :: dt_avg = 0

  real(r8), parameter :: cgs2mks = 1.0e-3_r8

  ! PIO descriptors (for restarts)

  type(var_desc_t), allocatable :: abstot_desc(:)
  type(var_desc_t) :: emstot_desc, absnxt_desc(4)

  logical  :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the zenith calculation
  real(r8) :: rad_uniform_angle = -99._r8

contains

  subroutine radiation_readnl(nlfile)

    use namelist_utils, only: find_group_name
    use units         , only: getunit, freeunit
    use spmd_utils    , only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, &
                               mpi_character, mpi_real8

    character(*), intent(in) :: nlfile

    integer unitn, ierr
    integer dtime      ! timestep size
    character(*), parameter :: sub = 'radiation_readnl'

    namelist /radiation_nl/ absems_data, iradsw, iradlw, iradae, irad_always, &
                            use_rad_dt_cosz, use_rad_uniform_angle, rad_uniform_angle

    if (masterproc) then
      unitn = getunit()
      open(unitn, file=trim(nlfile), status='old')
      call find_group_name(unitn, 'radiation_nl', status=ierr)
      if (ierr == 0) then
        read(unitn, radiation_nl, iostat=ierr)
        if (ierr /= 0) then
          call endrun(sub // ':: ERROR reading namelist')
        end if
      end if
      close(unitn)
      call freeunit(unitn)
    end if

    call mpi_bcast(absems_data          , len(absems_data), mpi_character, mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: absems_data')
    call mpi_bcast(iradsw               , 1               , mpi_integer  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: iradsw')
    call mpi_bcast(iradlw               , 1               , mpi_integer  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: iradlw')
    call mpi_bcast(iradae               , 1               , mpi_integer  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: iradae')
    call mpi_bcast(irad_always          , 1               , mpi_integer  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: irad_always')
    call mpi_bcast(use_rad_dt_cosz      , 1               , mpi_logical  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: use_rad_dt_cosz')
    call mpi_bcast(use_rad_uniform_angle, 1               , mpi_logical  , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: use_rad_uniform_angle')
    call mpi_bcast(rad_uniform_angle    , 1               , mpi_real8    , mstrid, mpicom, ierr)
    if (ierr /= 0) call endrun(sub // ': FATAL: mpi_bcast: rad_uniform_angle')

    if (use_rad_uniform_angle .and. rad_uniform_angle == -99.0_r8) then
      call endrun(sub // ' ERROR - use_rad_uniform_angle is set to .true, but rad_uniform_angle is not set!')
    end if

    ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary
    dtime  = get_step_size()
    if (iradsw      < 0) iradsw      = nint((-iradsw      * 3600) / dtime)
    if (iradlw      < 0) iradlw      = nint((-iradlw      * 3600) / dtime)
    if (irad_always < 0) irad_always = nint((-irad_always * 3600) / dtime)

    ! Convert iradae from hours to timesteps if necessary and check that
    ! iradae must be an even multiple of iradlw
    if (iradae < 0) iradae = nint((-iradae * 3600) / dtime)
    if (mod(iradae, iradlw) /= 0) then
      write(iulog, *) sub//': iradae must be an even multiple of iradlw.'
      write(iulog, *) '     iradae = ',iradae,', iradlw = ',iradlw
      call endrun(sub // ': iradae must be an even multiple of iradlw.')
    end if

   !----------------------------------------------------------------------- 
   ! Print runtime options to log.
   !-----------------------------------------------------------------------

   if (masterproc) then
      write(iulog,  *) 'CAMRT radiation scheme parameters:'
      write(iulog, 10) iradsw, iradlw, iradae, irad_always, use_rad_dt_cosz
      write(iulog,  *) '  Abs/Emis dataset: ', trim(absems_data)
   end if

10 format('  Frequency (timesteps) of Shortwave Radiation calc:    ',i5/, &
          '  Frequency (timesteps) of Longwave Radiation calc:     ',i5/, &
          '  Frequency (timesteps) of Absorptivity/Emissivity calc:',i5/, &
          '  SW/LW calc done every timestep for first N steps. N=  ',i5/, &
          '  Use average zenith angle:                             ',l5)

  end subroutine radiation_readnl

  subroutine radiation_register

    use physics_buffer, only: pbuf_add_field, dtype_r8
    use radiation_data, only: rad_data_register

    call pbuf_add_field('QRS' , 'global', dtype_r8, [pcols,pver], qrs_idx ) ! Shortwave radiative heating rate 
    call pbuf_add_field('QRL' , 'global', dtype_r8, [pcols,pver], qrl_idx ) ! Longwave  radiative heating rate 
    call pbuf_add_field('FSDS', 'global', dtype_r8, [pcols     ], fsds_idx) ! Surface solar downward flux
    call pbuf_add_field('FSNS', 'global', dtype_r8, [pcols     ], fsns_idx) ! Surface net shortwave flux
    call pbuf_add_field('FSNT', 'global', dtype_r8, [pcols     ], fsnt_idx) ! Top-of-model net shortwave flux
    call pbuf_add_field('FLNS', 'global', dtype_r8, [pcols     ], flns_idx) ! Surface net longwave flux
    call pbuf_add_field('FLNT', 'global', dtype_r8, [pcols     ], flnt_idx) ! Top-of-model net longwave flux

    call rad_data_register()

  end subroutine radiation_register

  pure logical function radiation_do(op, timestep) result(res)

    character(*), intent(in) :: op ! Name of operation
    integer, intent(in), optional :: timestep

    integer nstep ! Current timestep number

    if (present(timestep)) then
      nstep = timestep
    else
      nstep = get_nstep()
    end if

    select case (op)
    case ('sw') ! do a shortwave heating calc this timestep?
      res = nstep == 0 .or. iradsw == 1                     &
          .or. (mod(nstep-1, iradsw) == 0 .and. nstep /= 1) &
          .or. nstep <= irad_always
    case ('lw') ! do a longwave heating calc this timestep?
      res = nstep == 0 .or. iradlw == 1                     &
          .or. (mod(nstep-1, iradlw) == 0 .and. nstep /= 1) &
          .or. nstep <= irad_always
    case ('absems') ! do an absorptivity/emissivity calculation this timestep?
      res = nstep == 0 .or. iradae == 1                     &
          .or. (mod(nstep-1, iradae) == 0 .and. nstep /= 1)
    case ('aeres') ! write absorptivity/emissivity to restart file this timestep?
      res = mod(nstep, iradae) /= 0
    case default
      call endrun('radiation_do: unknown operation: ' // op)
    end select

  end function radiation_do

  real(r8) function radiation_nextsw_cday()
  
    ! Returns calendar day of next sw radiation calculation

    integer  nstep      ! timestep counter
    logical  dosw       ! true => do shosrtwave calc   
    integer  offset     ! offset for calendar day calculation
    integer  dTime      ! integer timestep size
    real(r8) calday     ! calendar day of 

    radiation_nextsw_cday = -1
    dosw   = .false.
    nstep  = get_nstep()
    dtime  = get_step_size()
    offset = 0
    do while (.not. dosw)
      nstep = nstep + 1
      offset = offset + dtime
      if (radiation_do('sw', nstep)) then
        radiation_nextsw_cday = get_curr_calday(offset=offset) 
        dosw = .true.
      end if
    end do
    if (radiation_nextsw_cday == -1) then
      call endrun('error in radiation_nextsw_cday')
    end if

  end function radiation_nextsw_cday

  subroutine radiation_init(pbuf2d)

    ! Initialize the radiation parameterization, add fields to the history buffer

    use cam_history   , only: addfld, add_default, horiz_only
    use physconst     , only: gravit, cpair, epsilo, stebol, &
                              pstd, mwdry, mwco2, mwo3
    use physics_buffer, only: physics_buffer_desc, pbuf_get_index
    use radsw         , only: radsw_init
    use radlw         , only: radlw_init
    use radae         , only: radae_init
    use radconstants  , only: radconstants_init
    use rad_solar_var , only: rad_solar_var_init
    use radiation_data, only: rad_data_init
    use phys_control  , only: phys_getopts
    use time_manager  , only: get_step_size

    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    integer nstep                       ! current timestep number
    logical history_amwg                ! output the variables used by the AMWG diag package
    logical history_vdiag               ! output the variables used by the AMWG variability diag package
    logical history_budget              ! output tendencies and state variables for CAM4
                                        ! temperature, water vapor, cloud ice and cloud
                                        ! liquid budgets.
    integer history_budget_histfile_num ! output history file number for budget fields
    integer dtime

    call radconstants_init()
    call rad_solar_var_init()

    call radsw_init(gravit)
    call radlw_init(gravit, stebol)
    call radae_init(gravit, epsilo, stebol, pstd, mwdry, mwco2, mwo3, absems_data)

    call rad_data_init(pbuf2d)

    ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation
    if (use_rad_dt_cosz)  then
      dtime  = get_step_size()
      dt_avg = real(iradsw * dtime, r8)
    end if

    ! Get physics buffer indices
    cld_idx = pbuf_get_index('CLD')
    rel_idx = pbuf_get_index('REL')
    rei_idx = pbuf_get_index('REI')

    ! "irad_always" is number of time steps to execute radiation continuously from start of
    ! initial OR restart run
    nstep = get_nstep()
    if (irad_always > 0) then
      nstep       = get_nstep()
      irad_always = irad_always + nstep
    end if

    ! Shortwave radiation
    call addfld('SOLIN'   , horiz_only, 'A', 'W/m2', 'Solar insolation'                               , sampling_seq='rad_lwsw')
    call addfld('SOLL'    , horiz_only, 'A', 'W/m2', 'Solar downward near infrared direct to surface' , sampling_seq='rad_lwsw')
    call addfld('SOLS'    , horiz_only, 'A', 'W/m2', 'Solar downward visible direct to surface'       , sampling_seq='rad_lwsw')
    call addfld('SOLLD'   , horiz_only, 'A', 'W/m2', 'Solar downward near infrared diffuse to surface', sampling_seq='rad_lwsw')
    call addfld('SOLSD'   , horiz_only, 'A', 'W/m2', 'Solar downward visible diffuse to surface'      , sampling_seq='rad_lwsw')
    call addfld('QRS'     , [ 'lev' ] , 'A', 'K/s' , 'Solar heating rate'                             , sampling_seq='rad_lwsw')
    call addfld('QRSC'    , [ 'lev' ] , 'A', 'K/s' , 'Clearsky solar heating rate'                    , sampling_seq='rad_lwsw')
    call addfld('FSNS'    , horiz_only, 'A', 'W/m2', 'Net solar flux at surface'                      , sampling_seq='rad_lwsw')
    call addfld('FSNT'    , horiz_only, 'A', 'W/m2', 'Net solar flux at top of model'                 , sampling_seq='rad_lwsw')
    call addfld('FSNTOA'  , horiz_only, 'A', 'W/m2', 'Net solar flux at top of atmosphere'            , sampling_seq='rad_lwsw')
    call addfld('FSUTOA'  , horiz_only, 'A', 'W/m2', 'Upwelling solar flux at top of atmosphere'      , sampling_seq='rad_lwsw')
    call addfld('FSNTOAC' , horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of atmosphere'   , sampling_seq='rad_lwsw')
    call addfld('FSDTOA'  , horiz_only, 'A', 'W/m2', 'Downwelling solar flux at top of atmosphere'    , sampling_seq='rad_lwsw')
    call addfld('FSN200'  , horiz_only, 'A', 'W/m2', 'Net shortwave flux at 200 mb'                   , sampling_seq='rad_lwsw')
    call addfld('FSN200C' , horiz_only, 'A', 'W/m2', 'Clearsky net shortwave flux at 200 mb'          , sampling_seq='rad_lwsw')
    call addfld('FSNTC'   , horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of model'        , sampling_seq='rad_lwsw')
    call addfld('FSNSC'   , horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at surface'             , sampling_seq='rad_lwsw')
    call addfld('FSDSC'   , horiz_only, 'A', 'W/m2', 'Clearsky downwelling solar flux at surface'     , sampling_seq='rad_lwsw')
    call addfld('FSDS'    , horiz_only, 'A', 'W/m2', 'Downwelling solar flux at surface'              , sampling_seq='rad_lwsw')
    call addfld('FUS'     , [ 'ilev' ], 'I', 'W/m2', 'Shortwave upward flux'                                                   )
    call addfld('FDS'     , [ 'ilev' ], 'I', 'W/m2', 'Shortwave downward flux'                                                 )
    call addfld('FUSC'    , [ 'ilev' ], 'I', 'W/m2', 'Shortwave clear-sky upward flux'                                         )
    call addfld('FDSC'    , [ 'ilev' ], 'I', 'W/m2', 'Shortwave clear-sky downward flux'                                       )
    call addfld('FSNIRTOA', horiz_only, 'A', 'W/m2', 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere'         , sampling_seq='rad_lwsw')
    call addfld('FSNRTOAC', horiz_only, 'A', 'W/m2', 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw')
    call addfld('FSNRTOAS', horiz_only, 'A', 'W/m2', 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere'        , sampling_seq='rad_lwsw')
    call addfld('FSNR'    , horiz_only, 'A', 'W/m2', 'Net solar flux at tropopause'                                        , sampling_seq='rad_lwsw')
    call addfld('SWCF'    , horiz_only, 'A', 'W/m2', 'Shortwave cloud forcing'                                             , sampling_seq='rad_lwsw')

    call addfld('TOT_CLD_VISTAU' , [ 'lev' ], 'A', '1', 'Total gbx cloud visible sw optical depth', sampling_seq='rad_lwsw', flag_xyfill=.true.)
    call addfld('TOT_ICLD_VISTAU', [ 'lev' ], 'A', '1', 'Total in-cloud visible sw optical depth' , sampling_seq='rad_lwsw', flag_xyfill=.true.)
    call addfld('LIQ_ICLD_VISTAU', [ 'lev' ], 'A', '1', 'Liquid in-cloud visible sw optical depth', sampling_seq='rad_lwsw', flag_xyfill=.true.)
    call addfld('ICE_ICLD_VISTAU', [ 'lev' ], 'A', '1', 'Ice in-cloud visible sw optical depth'   , sampling_seq='rad_lwsw', flag_xyfill=.true.)

    ! Longwave radiation
    call addfld('QRL'            , [ 'lev' ] , 'A', 'K/s' , 'Longwave heating rate'                           , sampling_seq='rad_lwsw')
    call addfld('QRLC'           , [ 'lev' ] , 'A', 'K/s' , 'Clearsky longwave heating rate'                  , sampling_seq='rad_lwsw')
    call addfld('FLNS'           , horiz_only, 'A', 'W/m2', 'Net longwave flux at surface'                    , sampling_seq='rad_lwsw')
    call addfld('FLDS'           , horiz_only, 'A', 'W/m2', 'Downwelling longwave flux at surface'            , sampling_seq='rad_lwsw')
    call addfld('FLNT'           , horiz_only, 'A', 'W/m2', 'Net longwave flux at top of model'               , sampling_seq='rad_lwsw')
    call addfld('FLUT'           , horiz_only, 'A', 'W/m2', 'Upwelling longwave flux at top of model'         , sampling_seq='rad_lwsw')
    call addfld('FLUTC'          , horiz_only, 'A', 'W/m2', 'Clearsky upwelling longwave flux at top of model', sampling_seq='rad_lwsw')
    call addfld('FLNTC'          , horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at top of model'      , sampling_seq='rad_lwsw')
    call addfld('FLN200'         , horiz_only, 'A', 'W/m2', 'Net longwave flux at 200 mb'                     , sampling_seq='rad_lwsw')
    call addfld('FLN200C'        , horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at 200 mb'            , sampling_seq='rad_lwsw')
    call addfld('FLNR'           , horiz_only, 'A', 'W/m2', 'Net longwave flux at tropopause'                 , sampling_seq='rad_lwsw')
    call addfld('FLNSC'          , horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at surface'           , sampling_seq='rad_lwsw')
    call addfld('FLDSC'          , horiz_only, 'A', 'W/m2', 'Clearsky downwelling longwave flux at surface'   , sampling_seq='rad_lwsw')
    call addfld('LWCF'           , horiz_only, 'A', 'W/m2', 'Longwave cloud forcing'                          , sampling_seq='rad_lwsw')
    call addfld('FUL'            , [ 'ilev' ], 'I', 'W/m2', 'Longwave upward flux'                                                     )
    call addfld('FDL'            , [ 'ilev' ], 'I', 'W/m2', 'Longwave downward flux'                                                   )
    call addfld('FULC'           , [ 'ilev' ], 'I', 'W/m2', 'Longwave clear-sky upward flux'                                           )
    call addfld('FDLC'           , [ 'ilev' ], 'I', 'W/m2', 'Longwave clear-sky downward flux'                                         )

    ! Heating rate needed for d(theta)/dt computation
    call addfld('HR',            , [ 'lev' ],  'A','K/s', 'Heating rate needed for d(theta)/dt computation')
   
    ! determine default variables
    call phys_getopts(history_amwg_out                =history_amwg,   &
                      history_vdiag_out               =history_vdiag,  &
                      history_budget_out              =history_budget, &
                      history_budget_histfile_num_out =history_budget_histfile_num)

    if (history_amwg) then
      ! Shortwave variables
      call add_default('SOLIN   ', 1, ' ')
      call add_default('QRS     ', 1, ' ')
      call add_default('FSNS    ', 1, ' ')
      call add_default('FSNT    ', 1, ' ')
      call add_default('FSDTOA  ', 1, ' ')
      call add_default('FSNTOA  ', 1, ' ')
      call add_default('FSUTOA  ', 1, ' ')
      call add_default('FSNTOAC ', 1, ' ')
      call add_default('FSNTC   ', 1, ' ')
      call add_default('FSNSC   ', 1, ' ')
      call add_default('FSDSC   ', 1, ' ')
      call add_default('FSDS    ', 1, ' ')
      call add_default('SWCF    ', 1, ' ')
      ! Longwave variables
      call add_default('QRL     ', 1, ' ')
      call add_default('FLNS    ', 1, ' ')
      call add_default('FLDS    ', 1, ' ')
      call add_default('FLNT    ', 1, ' ')
      call add_default('FLUT    ', 1, ' ')
      call add_default('FLUTC   ', 1, ' ')
      call add_default('FLNTC   ', 1, ' ')
      call add_default('FLNSC   ', 1, ' ')
      call add_default('FLDSC   ', 1, ' ')
      call add_default('LWCF    ', 1, ' ')
    end if
    if (single_column .and. scm_crm_mode) then
      ! Shortwave variables
      call add_default('FUS     ', 1, ' ')
      call add_default('FUSC    ', 1, ' ')
      call add_default('FDS     ', 1, ' ')
      call add_default('FDSC    ', 1, ' ')
      ! Longwave variables
      call add_default('FUL     ', 1, ' ')
      call add_default('FULC    ', 1, ' ')
      call add_default('FDL     ', 1, ' ')
      call add_default('FDLC    ', 1, ' ')
    end if
    
    if (history_budget .and. history_budget_histfile_num > 1) then
      call add_default('QRL     ', history_budget_histfile_num, ' ')
      call add_default('QRS     ', history_budget_histfile_num, ' ')
    end if
 
    if (history_vdiag) then
      call add_default('FLUT', 2, ' ')
      call add_default('FLUT', 3, ' ')
    end if
   
    cicewp_idx  = pbuf_get_index('CICEWP' )
    cliqwp_idx  = pbuf_get_index('CLIQWP' )
    cldemis_idx = pbuf_get_index('CLDEMIS')
    cldtau_idx  = pbuf_get_index('CLDTAU' )
    nmxrgn_idx  = pbuf_get_index('NMXRGN' )
    pmxrgn_idx  = pbuf_get_index('PMXRGN' )

  end subroutine radiation_init

  subroutine radiation_define_restart(file)

    type(file_desc_t), intent(inout) :: file

    integer i, ierr
    integer grid_id
    integer hdimcnt
    integer pver_id, pverp_id
    integer vsize
    integer dimids(4)

    type(cam_grid_header_info_t) :: info

    character(16) pname

    call pio_seterrorhandling(File, PIO_BCAST_ERROR)

    if (radiation_do('aeres')) then
      grid_id = cam_grid_id('physgrid')
      call cam_grid_write_attr(File, grid_id, info)
      hdimcnt = info%num_hdims()
      do i = 1, hdimcnt
        dimids(i) = info%get_hdimid(i)
      end do

      call cam_pio_def_dim(File, 'lev',  pver,  pver_id,  existOK=.true.)
      call cam_pio_def_dim(File, 'ilev', pverp, pverp_id, existOK=.true.)

      vsize = pverp - ntoplw + 1
      if (vsize /= pverp) then
        ierr = pio_def_dim(File, 'lwcols', vsize, dimids(hdimcnt+1))
      else
        dimids(hdimcnt+1) = pverp_id
      end if

      ! Split into vsize variables to avoid excessive memory usage in IO

      allocate(abstot_desc(ntoplw:pverp))

      do i = ntoplw, pverp
        write(pname, '(a,i3.3)') 'NAL_absorp', i
        ierr = pio_def_var(File, trim(pname), pio_double, dimids(1:hdimcnt+1), abstot_desc(i))
      end do
	
      dimids(hdimcnt+1) = pverp_id
      ierr = pio_def_var(File, 'Emissivity', pio_double, dimids(1:hdimcnt+1), emstot_desc)

      dimids(hdimcnt+1) = pver_id
      do i = 1, 4
        write(pname, '(a,i3.3)') 'NN_absorp', i
        ierr = pio_def_var(File, pname, pio_double, dimids(1:hdimcnt+1), absnxt_desc(i))
      end do
    end if

  end subroutine radiation_define_restart
  
  subroutine radiation_write_restart(file)

    type(file_desc_t), intent(inout) :: file

    integer i, ierr
    integer physgrid
    integer dims(3), gdims(3)
    integer nhdims
    integer ncol

    if (radiation_do('aeres')) then
      physgrid = cam_grid_id('physgrid')
      call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)

      do i = begchunk, endchunk
        ncol = get_ncols_p(i)
        if (ncol < pcols) then
          abstot_3d(ncol+1:pcols,:,:,i) = fillvalue
          absnxt_3d(ncol+1:pcols,:,:,i) = fillvalue
          emstot_3d(ncol+1:pcols,:,i)   = fillvalue
        end if
      end do
      
      ! abstot_3d is written as a series of 3D variables

      dims(1) = size(abstot_3d, 1) ! Should be pcols
      dims(2) = size(abstot_3d, 2) ! Should be (pverp-ntoplw+1)
      dims(3) = size(abstot_3d, 4) ! Should be endchunk - begchunk + 1
      gdims(nhdims+1) = dims(2)
      do i = ntoplw, pverp
        call cam_grid_write_dist_array(File, physgrid, dims(1:3), &
          gdims(1:nhdims+1), abstot_3d(:,:,i,:), abstot_desc(i))
      end do

      dims(1) = size(emstot_3d, 1) ! Should be pcols
      dims(2) = size(emstot_3d, 2) ! Should be pverp
      dims(3) = size(emstot_3d, 3) ! Should be endchunk - begchunk + 1
      gdims(nhdims+1) = dims(2)
      call cam_grid_write_dist_array(File, physgrid, dims(1:3), &
        gdims(1:nhdims+1), emstot_3d, emstot_desc)

      dims(1) = size(absnxt_3d, 1) ! Should be pcols
      dims(2) = size(absnxt_3d, 2) ! Should be pver
      dims(3) = size(absnxt_3d, 4) ! Should be endchunk - begchunk + 1
      gdims(nhdims+1) = dims(2)
      do i = 1, 4
        call cam_grid_write_dist_array(File, physgrid, dims(1:3), &
          gdims(1:nhdims+1), absnxt_3d(:,:,i,:), absnxt_desc(i))
      end do

      ! module data was allocated in radiation_define_restart
      deallocate(abstot_desc)
    end if

  end subroutine radiation_write_restart
  
  subroutine radiation_read_restart(file)

    type(file_desc_t), intent(inout) :: file

    integer err_handling
    integer ierr
    integer physgrid
    integer dims(3), gdims(3), nhdims
    integer vsize
    integer i

    type(var_desc_t) vardesc
    character(16) pname

    ! Put this call here for now.  It should move to an init method when the
    ! initialization and restart sequencing is unified.
    call initialize_radbuffer()

    if (radiation_do('aeres')) then
      call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling)
      ierr = pio_inq_varid(File, 'Emissivity', vardesc)
      call pio_seterrorhandling(File, err_handling)
      if (ierr /= PIO_NOERR) then
        if (masterproc) write(iulog,*) 'Warning: Emissivity variable not found on restart file.'
        return
      end if

      physgrid = cam_grid_id('physgrid')
      call cam_grid_dimensions(physgrid, gdims(1:2), nhdims)

      dims(1) = pcols
      dims(2) = pverp
      dims(3) = endchunk - begchunk + 1
      gdims(nhdims+1) = dims(2)

      call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
        gdims(1:nhdims+1), emstot_3d, vardesc)
        
      vsize = pverp - ntoplw + 1
      dims(2) = vsize
      gdims(nhdims+1) = dims(2)
        
      do i = ntoplw, pverp
        write(pname, '(A, I3.3)') 'NAL_absorp', i
        ierr = pio_inq_varid(File, trim(pname), vardesc)
        call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
          gdims(1:nhdims+1), abstot_3d(:,:,i,:), vardesc)
      end do

      dims(2) = pver
      gdims(nhdims+1) = dims(2)
      do i = 1, 4
        write(pname,'(a,i3.3)') 'NN_absorp', i
        ierr = pio_inq_varid(File, trim(pname), vardesc)
        call cam_grid_read_dist_array(File, physgrid, dims(1:3), &
          gdims(1:nhdims+1), absnxt_3d(:,:,i,:), vardesc)
      end do
    end if

  end subroutine radiation_read_restart
  
  subroutine radiation_tend(state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out)

    ! NOTE: Radiation uses cgs units, so conversions must be done from
    !       model fields to radiation fields.

    use physics_buffer   , only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx
    use phys_grid        , only: get_rlat_all_p, get_rlon_all_p
    use physics_types    , only: physics_state, physics_ptend
    use time_manager     , only: get_curr_calday
    use radheat          , only: radheat_tend
    use physconst        , only: cpair, stebol
    use radconstants     , only: nlwbands, nswbands
    use radsw            , only: radcswmx
    use radlw            , only: radclwmx
    use rad_constituents , only: rad_cnst_get_gas, rad_cnst_out
    use aer_rad_props    , only: aer_rad_props_sw, aer_rad_props_lw
    use interpolate_data , only: vertinterp
    use radiation_data   , only: rad_data_write
    use cloud_cover_diags, only: cloud_cover_diags_out
    use tropopause       , only: tropopause_find, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
    use orbit            , only: zenith

    type(physics_state),       target, intent(in)    :: state
    type(physics_ptend),               intent(out)   :: ptend
    type(physics_buffer_desc),         pointer       :: pbuf(:)
    type(cam_out_t),                   intent(inout) :: cam_out
    type(cam_in_t),                    intent(in)    :: cam_in
    real(r8),                          intent(out)   :: net_flx(pcols)
    type(rad_out_t), target, optional, intent(out)   :: rd_out

    type(rad_out_t), pointer :: rd  ! Allow rd_out to be optional by allocating a local object
                                    ! if the argument is not present

    integer i, k, lchnk, ncol
    logical dosw, dolw, doabsems
    integer , pointer :: nmxrgn(:)              ! pbuf pointer to Number of maximally overlapped regions
    real(r8), pointer :: pmxrgn(:,:)            ! Maximum values of pressure for each
                                                !    maximally overlapped region.
                                                !    0->pmxrgn(i,1) is range of pressure for
                                                !    1st region,pmxrgn(i,1)->pmxrgn(i,2) for
                                                !    2nd region, etc

    real(r8), pointer :: emis  (:,:)            ! Cloud longwave emissivity
    real(r8), pointer :: cldtau(:,:)            ! Cloud longwave optical depth
    real(r8), pointer :: cicewp(:,:)            ! in-cloud cloud ice water path
    real(r8), pointer :: cliqwp(:,:)            ! in-cloud cloud liquid water path

    real(r8) cltot(pcols)                       ! Diagnostic total cloud cover
    real(r8) cllow(pcols)                       !       "     low  cloud cover
    real(r8) clmed(pcols)                       !       "     mid  cloud cover
    real(r8) clhgh(pcols)                       !       "     hgh  cloud cover
    real(r8) ftem (pcols,pver)                  ! Temporary workspace for outfld variables

    integer itim_old
    real(r8), pointer :: rel(:,:)               ! Liquid effective drop radius (microns)
    real(r8), pointer :: rei(:,:)               ! Ice effective drop size (microns)
    real(r8), pointer :: cld(:,:)               ! Cloud fraction
    real(r8), pointer :: qrs(:,:)               ! Shortwave radiative heating rate
    real(r8), pointer :: qrl(:,:)               ! Longwave  radiative heating rate

    real(r8) calday                             ! Current calendar day
    real(r8) clat  (pcols)                      ! Current latitudes(radians)
    real(r8) clon  (pcols)                      ! Current longitudes(radians)
    real(r8) coszrs(pcols)                      ! Cosine solar zenith angle

    real(r8) fns   (pcols,pverp)                ! Net shortwave flux
    real(r8) fcns  (pcols,pverp)                ! Net clear-sky shortwave flux
    real(r8) fnl   (pcols,pverp)                ! Net longwave flux
    real(r8) fcnl  (pcols,pverp)                ! Net clear-sky longwave flux

    ! This is used by the chemistry.
    real(r8), pointer :: fsds(:)                ! Surface solar down flux

    ! This is used for the energy checker and the Eulerian dycore.
    real(r8), pointer :: fsns(:)                ! Surface solar absorbed flux
    real(r8), pointer :: fsnt(:)                ! Net column abs solar flux at model top
    real(r8), pointer :: flns(:)                ! Srf longwave cooling (up-down) flux
    real(r8), pointer :: flnt(:)                ! Net outgoing lw flux at model top

    real(r8) pbr    (pcols,pver)                ! Model mid-level pressures (dynes/cm2)
    real(r8) pnm    (pcols,pverp)               ! Model interface pressures (dynes/cm2)
    real(r8) eccf                               ! Earth/sun distance factor
    real(r8) lwupcgs(pcols)                     ! Upward longwave flux in CGS units
 
    real(r8), pointer :: n2o   (:,:)            ! Nitrous oxide mass mixing ratio
    real(r8), pointer :: ch4   (:,:)            ! Methane mass mixing ratio
    real(r8), pointer :: cfc11 (:,:)            ! CFC11 mass mixing ratio
    real(r8), pointer :: cfc12 (:,:)            ! CFC12 mass mixing ratio
    real(r8), pointer :: o3    (:,:)            ! Ozone mass mixing ratio
    real(r8), pointer :: o2    (:,:)            ! Oxygen mass mixing ratio
    real(r8)             o2_col(pcols)          ! Column oxygen mmr
    real(r8), pointer :: co2   (:,:)            ! CO2 mass mixing ratio
    real(r8)       co2_col_mean(pcols)          ! CO2 column mean mmr
    real(r8), pointer :: sp_hum(:,:)            ! Specific humidity

    ! Aerosol shortwave radiative properties
    real(r8) aer_tau    (pcols,0:pver,nswbands) ! Aerosol extinction optical depth
    real(r8) aer_tau_w  (pcols,0:pver,nswbands) ! Aerosol single scattering albedo * tau
    real(r8) aer_tau_w_g(pcols,0:pver,nswbands) ! Aerosol assymetry parameter * w * tau
    real(r8) aer_tau_w_f(pcols,0:pver,nswbands) ! Aerosol forward scattered fraction * w * tau

    ! Aerosol longwave absorption optical depth
    real(r8) odap_aer(pcols,pver,nlwbands)

    ! Gathered indicies of day and night columns
    !  chunk_column_index = IdxDay(daylight_column_index)
    integer Nday                                ! Number of daylight columns
    integer Nnite                               ! Number of night columns
    integer :: IdxDay (pcols)                   ! Indicies of daylight coumns
    integer :: IdxNite(pcols)                   ! Indicies of night coumns

    character(*), parameter :: name = 'radiation_tend'

    ! Tropopause diagnostic
    integer  troplev(pcols)
    real(r8) p_trop (pcols)

    logical write_output ! switch for outfld calls

    lchnk = state%lchnk
    ncol  = state%ncol

    calday = get_curr_calday()

    if (present(rd_out)) then
      rd => rd_out
      write_output = .false.
    else
      allocate(rd)
      write_output=.true.
    end if

    itim_old = pbuf_old_tim_idx()
    call pbuf_get_field(pbuf, cld_idx , cld , start=[1,1,itim_old], kount=[pcols,pver,1])
    call pbuf_get_field(pbuf, qrs_idx , qrs )
    call pbuf_get_field(pbuf, qrl_idx , qrl )
    call pbuf_get_field(pbuf, fsds_idx, fsds)
    call pbuf_get_field(pbuf, fsns_idx, fsns)
    call pbuf_get_field(pbuf, fsnt_idx, fsnt)
    call pbuf_get_field(pbuf, flns_idx, flns)
    call pbuf_get_field(pbuf, flnt_idx, flnt)
    call pbuf_get_field(pbuf, rel_idx , rel )
    call pbuf_get_field(pbuf, rei_idx , rei )
   
    ! For CRM, make cloud equal to input observations:
    if (single_column .and. scm_crm_mode .and. have_cld) then
      do k = 1, pver
        cld(:ncol,k) = cldobs(k)
      end do
    end if

    ! Cosine solar zenith angle for current time step
    call get_rlat_all_p(lchnk, ncol, clat)
    call get_rlon_all_p(lchnk, ncol, clon)
    if (use_rad_uniform_angle) then
      call zenith(calday, clat, clon, coszrs, ncol, dt_avg, uniform_angle=rad_uniform_angle)
    else
      call zenith(calday, clat, clon, coszrs, ncol, dt_avg)
    end if

    ! Gather night/day column indices.
    Nday  = 0
    Nnite = 0
    do i = 1, ncol
      if (coszrs(i) > 0) then
        Nday = Nday + 1
        IdxDay(Nday) = i
      else
        Nnite = Nnite + 1
        IdxNite(Nnite) = i
      end if
    end do

    dosw     = radiation_do('sw')      ! do shortwave heating calc this timestep?
    dolw     = radiation_do('lw')      ! do longwave heating calc this timestep?
    doabsems = radiation_do('absems')  ! do absorptivity/emissivity calc this timestep?

    if (dosw .or. dolw) then
      ! pbuf cloud properties set in cloud_diagnostics
      call pbuf_get_field(pbuf, cicewp_idx , cicewp)
      call pbuf_get_field(pbuf, cliqwp_idx , cliqwp)
      call pbuf_get_field(pbuf, cldemis_idx, emis  )
      call pbuf_get_field(pbuf, cldtau_idx , cldtau)
      call pbuf_get_field(pbuf, pmxrgn_idx , pmxrgn)
      call pbuf_get_field(pbuf, nmxrgn_idx , nmxrgn)

      ! For CRM, make cloud liquid water path equal to input observations
      if (single_column .and. scm_crm_mode .and. have_clwp) then
        do k = 1, pver
          cliqwp(:ncol,k) = clwpobs(k)
        end do
      end if

      ! Get specific humidity
      call rad_cnst_get_gas(0, 'H2O', state, pbuf,  sp_hum)

      ! Get ozone mass mixing ratio.
      call rad_cnst_get_gas(0, 'O3' , state, pbuf,  o3)

      ! Get CO2 mass mixing ratio and compute column mean values
      call rad_cnst_get_gas(0, 'CO2', state, pbuf,  co2)
      call calc_col_mean(state, co2, co2_col_mean)

      ! Construct CGS unit reps of pmid and pint and get "eccf" - earthsundistancefactor
      call radinp(ncol, state%pmid, state%pint, pbr, pnm, eccf)

      ! Solar radiation computation
      if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then
        call tropopause_find(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
      end if

      if (dosw) then
        ! Get Oxygen mass mixing ratio.
        call rad_cnst_get_gas(0, 'O2', state, pbuf,  o2)
        call calc_col_mean(state, o2, o2_col)
   
        ! Get aerosol radiative properties.
        call aer_rad_props_sw(0, state, pbuf,  nnite, idxnite, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f)

        call radcswmx(                      &
          lchnk        =lchnk             , &
          ncol         =ncol              , &
          E_pint       =pnm               , &
          E_pmid       =pbr               , &
          E_h2ommr     =sp_hum            , &
          E_o3mmr      =o3                , &
          E_o2mmr      =o2_col            , &
          E_cld        =cld               , &
          E_cicewp     =cicewp            , &
          E_cliqwp     =cliqwp            , &
          E_rel        =rel               , &
          E_rei        =rei               , &
          eccf         =eccf              , &
          E_coszrs     =coszrs            , &
          solin        =rd%solin          , &
          E_asdir      =cam_in%asdir      , &
          E_asdif      =cam_in%asdif      , &
          E_aldir      =cam_in%aldir      , &
          E_aldif      =cam_in%aldif      , &
          nmxrgn       =nmxrgn            , &
          pmxrgn       =pmxrgn            , &
          qrs          =qrs               , &
          qrsc         =rd%qrsc           , &
          fsnt         =fsnt              , &
          fsntc        =rd%fsntc          , &
          fsdtoa       =rd%fsdtoa         , &
          fsntoa       =rd%fsntoa         , &
          fsutoa       =rd%fsutoa         , &
          fsntoac      =rd%fsntoac        , &
          fsnirtoa     =rd%fsnirt         , &
          fsnrtoac     =rd%fsnrtc         , &
          fsnrtoaq     =rd%fsnirtsq       , &
          fsns         =fsns              , &
          fsnsc        =rd%fsnsc          , &
          fsdsc        =rd%fsdsc          , &
          fsds         =fsds              , &
          sols         =cam_out%sols      , &
          soll         =cam_out%soll      , &
          solsd        =cam_out%solsd     , &
          solld        =cam_out%solld     , &
          fns          =fns               , &
          fcns         =fcns              , &
          Nday         =Nday              , &
          Nnite        =Nnite             , &
          IdxDay       =IdxDay            , &
          IdxNite      =IdxNite           , &
          E_co2mmr     =co2_col_mean      , &
          E_aer_tau    =aer_tau           , &
          E_aer_tau_w  =aer_tau_w         , &
          E_aer_tau_w_g=aer_tau_w_g       , &
          E_aer_tau_w_f=aer_tau_w_f       , &
          tauxcl_out   =rd%liq_icld_vistau, &
          tauxci_out   =rd%ice_icld_vistau  ) 

        ! Output net fluxes at 200 mb
        call vertinterp(ncol, pcols, pverp, state%pint, 20000.0_r8, fcns, rd%fsn200c)
        call vertinterp(ncol, pcols, pverp, state%pint, 20000.0_r8, fns,  rd%fsn200)
        if (hist_fld_active('FSNR')) then
          do i = 1, ncol
            call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i))
            rd%fsnr(i) = rd%fsnr(i)*cgs2mks
          end do
        else
          rd%fsnr(:) = 0
        end if

        ! Convert units of shortwave fields needed by rest of model from CGS to MKS
        do i = 1, ncol
          rd%solin   (i) = rd%solin   (i) * cgs2mks
          fsds       (i) = fsds       (i) * cgs2mks
          rd%fsnirt  (i) = rd%fsnirt  (i) * cgs2mks
          rd%fsnrtc  (i) = rd%fsnrtc  (i) * cgs2mks
          rd%fsnirtsq(i) = rd%fsnirtsq(i) * cgs2mks
          fsnt       (i) = fsnt       (i) * cgs2mks
          rd%fsdtoa  (i) = rd%fsdtoa  (i) * cgs2mks
          fsns       (i) = fsns       (i) * cgs2mks
          rd%fsntc   (i) = rd%fsntc   (i) * cgs2mks
          rd%fsnsc   (i) = rd%fsnsc   (i) * cgs2mks
          rd%fsdsc   (i) = rd%fsdsc   (i) * cgs2mks
          rd%fsntoa  (i) = rd%fsntoa  (i) * cgs2mks
          rd%fsutoa  (i) = rd%fsutoa  (i) * cgs2mks
          rd%fsntoac (i) = rd%fsntoac (i) * cgs2mks
          rd%fsn200  (i) = rd%fsn200  (i) * cgs2mks
          rd%fsn200c (i) = rd%fsn200c (i) * cgs2mks
          rd%swcf    (i) = rd%fsntoa  (i) - rd%fsntoac(i)
        end do

        ! Initialize tau_cld_vistau and tau_icld_vistau as fillvalue, they will stay fillvalue for night columns
        rd%tot_icld_vistau(1:pcols,1:pver) = fillvalue
        rd%tot_cld_vistau (1:pcols,1:pver) = fillvalue

        ! Only do calcs for tot_cld_vistau and tot_icld_vistau on daytime columns
        do i = 1, Nday
          ! Sum the water and ice optical depths to get total in-cloud cloud optical depth
          rd%tot_icld_vistau(IdxDay(i),1:pver) = rd%liq_icld_vistau(IdxDay(i),1:pver) + &
                                                 rd%ice_icld_vistau(IdxDay(i),1:pver)

          ! Sum wat and ice, multiply by cloud fraction to get grid-box value
          rd%tot_cld_vistau (IdxDay(i),1:pver) = (rd%liq_icld_vistau(IdxDay(i),1:pver) + &
                                                  rd%ice_icld_vistau(IdxDay(i),1:pver))*cld(IdxDay(i),1:pver)
        end do
        ! Add fillvalue for night columns
        do i = 1, Nnite
          rd%liq_icld_vistau(IdxNite(i),:) = fillvalue
          rd%ice_icld_vistau(IdxNite(i),:) = fillvalue
        end do
        if (write_output) call radiation_output_sw(state, rd, cam_out, fsns, fsnt, fsds, qrs)
      end if ! dosw

      ! Longwave radiation computation
      if (dolw) then
        ! Convert upward longwave flux units to CGS
        do i = 1, ncol
          lwupcgs(i) = cam_in%lwup(i) * 1000
          if (single_column .and. scm_crm_mode .and. have_tg) &
            lwupcgs(i) = 1000 * stebol * tground(1)**4
        end do

        ! Get gas phase constituents.
        call rad_cnst_get_gas(0, 'N2O'  , state, pbuf, n2o  )
        call rad_cnst_get_gas(0, 'CH4'  , state, pbuf, ch4  )
        call rad_cnst_get_gas(0, 'CFC11', state, pbuf, cfc11)
        call rad_cnst_get_gas(0, 'CFC12', state, pbuf, cfc12)

        ! absems requires lw absorption optical depth and transmission through aerosols
        if (doabsems) call aer_rad_props_lw(0, state, pbuf, odap_aer)

        call radclwmx(              &
          lchnk     =lchnk        , &
          ncol      =ncol         , &
          doabsems  =doabsems     , &
          lwupcgs   =lwupcgs      , &
          tnm       =state%t      , &
          qnm       =sp_hum       , &
          o3        =o3           , &
          pmid      =pbr          , &
          pint      =pnm          , &
          pmln      =state%lnpmid , &
          piln      =state%lnpint , &
          n2o       =n2o          , &
          ch4       =ch4          , &
          cfc11     =cfc11        , &
          cfc12     =cfc12        , &
          cld       =cld          , &
          emis      =emis         , &
          pmxrgn    =pmxrgn       , &
          nmxrgn    =nmxrgn       , &
          qrl       =qrl          , &
          qrlc      =rd%qrlc      , &
          flns      =flns         , &
          flnt      =flnt         , &
          flnsc     =rd%flnsc     , &
          flntc     =rd%flntc     , &
          flwds     =cam_out%flwds, &
          fldsc     =rd%fldsc     , &
          flut      =rd%flut      , &
          flutc     =rd%flutc     , &
          fnl       =fnl          , &
          fcnl      =fcnl         , &
          co2mmr    =co2_col_mean , &
          odap_aer  =odap_aer     )

        !  Output fluxes at 200 mb
        call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl,  rd%fln200)
        call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c)
        if (hist_fld_active('FLNR')) then
          do i = 1, ncol
            call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i))
          end do
        else
          rd%flnr(:) = 0
        end if

        ! Convert units of longwave fields needed by rest of model from CGS to MKS
        do i = 1, ncol
          flnt         (i) = flnt         (i) * cgs2mks
          rd%flut      (i) = rd%flut      (i) * cgs2mks
          rd%flutc     (i) = rd%flutc     (i) * cgs2mks
          rd%lwcf      (i) = rd%flutc     (i) - rd%flut(i)
          flns         (i) = flns         (i) * cgs2mks
          rd%fldsc     (i) = rd%fldsc     (i) * cgs2mks
          rd%flntc     (i) = rd%flntc     (i) * cgs2mks
          rd%fln200    (i) = rd%fln200    (i) * cgs2mks
          rd%fln200c   (i) = rd%fln200c   (i) * cgs2mks
          rd%flnsc     (i) = rd%flnsc     (i) * cgs2mks
          cam_out%flwds(i) = cam_out%flwds(i) * cgs2mks
          rd%flnr      (i) = rd%flnr      (i) * cgs2mks
        end do
        if (write_output) call radiation_output_lw(state, rd, cam_out, flns, flnt, qrl)
      end if ! dolw

      ! Output aerosol mmr
      if (write_output) call rad_cnst_out(0, state, pbuf)

      ! Cloud cover diagnostics
      ! radsw can change pmxrgn and nmxrgn so cldsav needs to follow radsw
      if (write_output) call cloud_cover_diags_out(lchnk, ncol, cld, state%pmid, nmxrgn, pmxrgn)
    else   !  if (dosw .or. dolw) then
      ! Convert radiative heating rates from Q*dp to Q for energy conservation
      do k = 1, pver
        do i = 1, ncol
          qrs(i,k) = qrs(i,k) / state%pdel(i,k)
          qrl(i,k) = qrl(i,k) / state%pdel(i,k)
        end do
      end do
    end if !  if (dosw .or. dolw) then

    ! Output radiation inputs and resulting heating rates
    call rad_data_write(pbuf, state, cam_in, coszrs)

    ! Compute net radiative heating tendency
    call radheat_tend(state, pbuf,  ptend, qrl, qrs, fsns, &
                      fsnt, flns, flnt, cam_in%asdir, net_flx)

    if (write_output) then
      ! Compute heating rate for dtheta/dt
      do k = 1, pver
        do i = 1, ncol
          ftem(i,k) = (qrs(i,k) + qrl(i,k)) / cpair * (1.e5_r8 / state%pmid(i,k))**cappa
        end do
      end do
      call outfld('HR      ', ftem, pcols, lchnk)
    end if

    ! Convert radiative heating rates to Q*dp for energy conservation
    do k =1 , pver
      do i = 1, ncol
        qrs(i,k) = qrs(i,k) * state%pdel(i,k)
        qrl(i,k) = qrl(i,k) * state%pdel(i,k)
      end do
    end do

    cam_out%netsw(:ncol) = fsns(:ncol)

    if (.not. present(rd_out)) then
      deallocate(rd)
    end if

  end subroutine radiation_tend

  subroutine radiation_output_sw(state, rd, cam_out, fsns, fsnt, fsds, qrs)

    ! Dump shortwave radiation information to history buffer (diagnostics)

    type(physics_state), intent(in) :: state
    type(rad_out_t),     intent(in) :: rd
    type(cam_out_t),     intent(in) :: cam_out
    real(r8),            intent(in) :: fsns(pcols) ! Surface solar absorbed flux
    real(r8),            intent(in) :: fsnt(pcols) ! Net column abs solar flux at model top
    real(r8),            intent(in) :: fsds(pcols) ! Surface solar down flux
    real(r8),            pointer    :: qrs(:,:)    ! shortwave radiative heating rate

    integer  lchnk, ncol
    real(r8) ftem(pcols,pver)

    lchnk = state%lchnk
    ncol  = state%ncol

    ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair
    call outfld('QRS     ', ftem         , pcols, lchnk)
    ftem(:ncol,:pver) = rd%qrsc(:ncol,:pver)/cpair
    call outfld('QRSC    ', ftem         , pcols, lchnk)
    call outfld('SOLIN   ', rd%solin     , pcols, lchnk)
    call outfld('FSDS    ', fsds         , pcols, lchnk)
    call outfld('FSNIRTOA', rd%fsnirt    , pcols, lchnk)
    call outfld('FSNRTOAC', rd%fsnrtc    , pcols, lchnk)
    call outfld('FSNRTOAS', rd%fsnirtsq  , pcols, lchnk)
    call outfld('FSNT    ', fsnt         , pcols, lchnk)
    call outfld('FSDTOA  ', rd%fsdtoa    , pcols, lchnk)
    call outfld('FSNS    ', fsns         , pcols, lchnk)
    call outfld('FSNTC   ', rd%fsntc     , pcols, lchnk)
    call outfld('FSNSC   ', rd%fsnsc     , pcols, lchnk)
    call outfld('FSDSC   ', rd%fsdsc     , pcols, lchnk)
    call outfld('FSNTOA  ', rd%fsntoa    , pcols, lchnk)
    call outfld('FSUTOA  ', rd%fsutoa    , pcols, lchnk)
    call outfld('FSNTOAC ', rd%fsntoac   , pcols, lchnk)
    call outfld('SOLS    ', cam_out%sols , pcols, lchnk)
    call outfld('SOLL    ', cam_out%soll , pcols, lchnk)
    call outfld('SOLSD   ', cam_out%solsd, pcols, lchnk)
    call outfld('SOLLD   ', cam_out%solld, pcols, lchnk)
    call outfld('FSN200  ', rd%fsn200    , pcols, lchnk)
    call outfld('FSN200C ', rd%fsn200c   , pcols, lchnk)
    call outfld('FSNR'    , rd%fsnr      , pcols, lchnk)
    call outfld('SWCF    ', rd%swcf      , pcols, lchnk)

    call outfld('TOT_CLD_VISTAU    ', rd%tot_cld_vistau  , pcols, lchnk)
    call outfld('TOT_ICLD_VISTAU   ', rd%tot_icld_vistau , pcols, lchnk)
    call outfld('LIQ_ICLD_VISTAU   ', rd%liq_icld_vistau , pcols, lchnk)
    call outfld('ICE_ICLD_VISTAU   ', rd%ice_icld_vistau , pcols, lchnk)
   
  end subroutine radiation_output_sw

  subroutine radiation_output_lw(state,  rd, cam_out, flns, flnt, qrl)

    ! Dump longwave radiation information to history tape buffer (diagnostics)

    type(physics_state), intent(in) :: state
    type(rad_out_t),     intent(in) :: rd
    type(cam_out_t),     intent(in) :: cam_out
    real(r8),            intent(in) :: flns(pcols) ! Srf longwave cooling (up-down) flux
    real(r8),            intent(in) :: flnt(pcols) ! Net outgoing lw flux at model top
    real(r8),            pointer    :: qrl(:,:)    ! longwave  radiative heating rate

    integer lchnk, ncol
    real(r8) ftem(pcols,pver)

    lchnk = state%lchnk
    ncol  = state%ncol

    call outfld('QRL     ', qrl(:ncol,:) / cpair    , ncol , lchnk)
    call outfld('QRLC    ', rd%qrlc(:ncol,:) / cpair, ncol , lchnk)
    call outfld('FLNT    ', flnt                    , pcols, lchnk)
    call outfld('FLUT    ', rd%flut                 , pcols, lchnk)
    call outfld('FLUTC   ', rd%flutc                , pcols, lchnk)
    call outfld('FLNTC   ', rd%flntc                , pcols, lchnk)
    call outfld('FLNS    ', flns                    , pcols, lchnk)
    call outfld('FLDS    ', cam_out%flwds           , pcols, lchnk)
    call outfld('FLNSC   ', rd%flnsc                , pcols, lchnk)
    call outfld('FLDSC   ', rd%fldsc                , pcols, lchnk)
    call outfld('LWCF    ', rd%lwcf                 , pcols, lchnk)
    call outfld('FLN200  ', rd%fln200               , pcols, lchnk)
    call outfld('FLN200C ', rd%fln200c              , pcols, lchnk)
    call outfld('FLNR '   , rd%flnr                 , pcols, lchnk)

  end subroutine radiation_output_lw

  subroutine radinp(ncol, pmid, pint, pmidrd, pintrd, eccf)

    use shr_orb_mod
    use time_manager, only: get_curr_calday

    integer , intent(in ) :: ncol                ! Number of atmospheric columns
    real(r8), intent(in ) :: pmid  (pcols,pver ) ! Pressure at model mid-levels (pascals)
    real(r8), intent(in ) :: pint  (pcols,pverp) ! Pressure at model interfaces (pascals)
    real(r8), intent(out) :: pmidrd(pcols,pver ) ! Pressure at mid-levels (dynes/cm*2)
    real(r8), intent(out) :: pintrd(pcols,pverp) ! Pressure at interfaces (dynes/cm*2)
    real(r8), intent(out) :: eccf                ! Earth-sun distance factor

    integer i, k
    real(r8) calday       ! Current calendar day
    real(r8) delta        ! Solar declination angle

    calday = get_curr_calday()
    call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, delta, eccf)

    ! Convert pressure from pascals to dynes/cm2
    do k = 1, pver
      do i = 1, ncol
        pmidrd(i,k) = pmid(i,k) * 10
        pintrd(i,k) = pint(i,k) * 10
      end do
    end do
    do i = 1, ncol
      pintrd(i,pverp) = pint(i,pverp) * 10
    end do
   
  end subroutine radinp

  subroutine calc_col_mean(state, mmr_pointer, mean_value)

    use cam_logfile,  only: iulog

    type(physics_state),        intent(in)  :: state
    real(r8), dimension(:,:),   pointer     :: mmr_pointer  ! mass mixing ratio (lev)
    real(r8), dimension(pcols), intent(out) :: mean_value   ! column mean mmr

    integer i, k, ncol
    real(r8) ptot(pcols)

    ncol       = state%ncol
    mean_value = 0
    ptot       = 0

    do k = 1, pver
      do i = 1, ncol
        mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k)
        ptot      (i) = ptot      (i) + state%pdeldry(i,k)
       end do
    end do
    do i = 1, ncol
      mean_value(i) = mean_value(i) / ptot(i)
    end do

  end subroutine calc_col_mean

end module radiation

